/******************************************************************************

                                QSHARE.C

  Title:        Shared utilities in physical space
  Author:       Hanna Szoke
  Date:         Dec. 22, 1993 
  Function:     I/O and analysis of (binary) data in real space.

void allocNpArrays() 	- solve		
void initNparrays() 	- solve
void updateKnownValues()- solve		
int  checkSolutions()	- solve
void writeSol() 	- solve
void analyzEpixFile()	- back, maketar, qshare
void readEpixFile() 	- convert, distance, maketar, qshare, regrid, variance
void writeEpixFile()	- back, convert, maketar, qshare, variance
real inv_mean_wt()	- qshare
void analyzeTar()	- qshare

*****************************************************************************/

#include "speden.h"	/* includes util.h, cellparams.h, symmetry.h, dims.h */
#include "speden_Np.h"

	/* declarations of speden_Np.h variables */

real	*nump ;                 /* electrons/voxel, one iteration */
real	*snump ;		/* electrons/voxel, all iterations */
real  *totp ;              	/* running total electrons/voxel   */
real  *knownp ;               /* electrons/voxel corresponding to fc file */
real  *target ;               /* electrons/voxel corresponding to targets */
real  *weightp ;              /* weights for Np-space cost function */
real  *extra ;        	/* extra work array for CS constraint */

static	void	analyzeTar(real *, real *, int) ;	
static	real	inv_mean_wt(real *) ;
static  char	solution_filename[MAXSTRING] ;

		/************************************************************
               	Allocate & initialize arrays in Np space needed by solve
		(excluding those needed only by conjugate gradient solver
		and those in Npextended).
		************************************************************/
void allocNpArrays() 
{
     int	NN = 0 ;
     int	k ;

     strcpy(message, "allocNpArrays") ;
     knownp = (real *) e_malloc(Nptotal*sizeof(real), message) ;

     NN = Ntargets*Nptotal ;

     if (NN > 0) {
          target = (real *) e_malloc(NN*sizeof(real), message) ;
          weightp= (real *) e_malloc(NN*sizeof(real), message) ;
     }
     for (k = 0; k < Nconstraints; k++) 
	  if (con_type[k] == CS) {
               extra = (real *) e_malloc(Nptotal*sizeof(real), message) ;
	       break ;
          }
     return ;
}

		/************************************************************
   		Decide whether size of solutions warrants a continuation of 
		the program.
		************************************************************/

#define EPS_GO_ON     1.e-3           /* Criterion for convergence */

int checkSolutions()
{
     int        n ;
     real     snumpmax ;
     real     numpmax ;

       /************************************************************
       Find largest solution in the cumulative array
       and in the current iteration.  
       ************************************************************/

     snumpmax = fabs (*snump) ;

     for (n = 1; n < Nptotal; n++)
          if ( fabs (*(snump+n)) >  snumpmax) {
               snumpmax = fabs (*(snump+n)) ;
          }

     numpmax = fabs (*nump) ;

     for (n = 1; n < Nptotal; n++)
          if ( fabs (*(nump+n)) >  numpmax) 
               numpmax = fabs (*(nump+n)) ;

       /************************************************************
       If it's too small, quit (unlikely in real cases).        
       ************************************************************/

     if (numpmax < snumpmax * EPS_GO_ON) {
          sprintf(message, 
		 "Solutions converged to within %g of maximum electrons/voxel", 
		 EPS_GO_ON) ;
          printTwice(message) ;
          return (STOP) ;
     }
     else
          return(GO_ON) ;
}

		/************************************************************
		Initialize Np space arrays, reading files where necessary.
		initMinMaxType() & alloc_ccg_mem() are called by Back as well.
		************************************************************/
void	initNparrays() 
{
     void	setup_csind() ;
     void	initMinMaxType() ;
     void	reviseMinMaxp() ;
     void	readEpixFile(char *, real *) ; 
     void	prepare_relative_index_list() ;	
     void	alloc_ccg_mem() ;		/* in ccg.c */
     int	hrInitialize() ;		/* in highres.c */
     void	hrSetup() ;			/* in highres.c */
     void	hrMaskOut() ;			/* in highres.c */

     real	volcoef ;	/*  V / (2pi*eta*dr^2)^(3/2)  */
     real	*big ;
     int	n ;
     int	k, NN, maxNhr ;

     for (n = 0; n < Nptotal; n++) 
	  *(knownp+n) = 0 ;

       /************************************************************
       Read model; in correction mode,                        
       revise the lower (and upper) bounds on the optimizer.
       Special treatment of "empty" is not really needed, since
       readEpixFile() deals with it...
       Note that hrInitialize(knownp) DISCARDS all the high points
       in knownp by masking them out!!!
       ************************************************************/

     if (strcmp (md_filename, "empty") != 0) {
          printTwice("Reading model electron map ...") ;
          readEpixFile(md_filename, knownp) ;
     }
     else
	  printTwice("Model electron map file is empty") ;

     setup_csind() ;

     if (HighRes)  {

       /************************************************************
        At this point, we re-allocate the maximum possible size for
        knownp and transfer Nptotal to this "big" array; then we find
        the actual req'd size (Npextended) and eventually reallocate
        and transfer data back to knownp.
	Maybe use memcpy(big, knownp, Nptotal*sizeof(real))?
       ************************************************************/

          maxNhr = hrInitialize(knownp) ;	
          sprintf(message, "initNparrays") ;

	  big = (real *) e_malloc((Nptotal + maxNhr)*sizeof(real), message);
	 
          for (n = 0; n < Nptotal; n++) 
	       *(big+n) = *(knownp+n) ;

          e_free(Nptotal*sizeof(real), knownp) ;

	  hrSetup(big) ;		

	  knownp = (real *) e_malloc(Npextended * sizeof(real), message);
	 
          for (n = 0; n < Npextended; n++) 
	       *(knownp+n) = *(big+n) ;

          e_free((Nptotal + maxNhr)*sizeof(real), big) ;

          hrMaskOut(knownp) ;
          sprintf(message, 
                "Nptotal = %d, Npextended = %d", Nptotal, Npextended) ;
          printTwice(message) ;
     }
     else{


     Npextended = Nptotal ;
 
     printTwice(
     "Mar 27, 2006: this code not yet corrected for complex density!/n") ;

     sprintf(message, "initNparrays") ;
     nump = (real *) e_malloc(Npextended*sizeof(real), message) ;
     snump  = (real *) e_malloc(Npextended*sizeof(real), message) ;
     totp   = (real *) e_malloc(Npextended*sizeof(real), message) ;

     for (n = 0; n < Npextended; n++) 
	  *(totp+n)  = *(knownp+n) ;
     
     initMinMaxType(Npextended) ;
     reviseMinMaxp(Npextended, knownp) ;

     alloc_ccg_mem(Npextended) ;

       /************************************************************
       Calculate a generic weight-based coefficient, volcoef. 
       ************************************************************/

     volcoef = crys_vol / pow((TWOPI*eta*dr*dr), (real) 1.5) ;
     sprintf(message, 
     "The volume coefficient: V/(2pi*eta*drsq)^3/2 = %g", volcoef) ;
     printTwice(message) ;

     for (k = 0, NN = 0; k < Nconstraints; k++)  {

       /************************************************************
       If there are targets, including phase extension targets,
       read them and their associated weights;  for ordinary
       targets only, factor the inverse-weight-based coefficient 
       into con_coef[k].
       Note: you may invoke "full" as a keyword to use default 
       weights (1 everywhere), rather than reading in a weight file.
       Report average el/voxel over target. 
       ************************************************************/

          if ((con_type[k] == TARGET)||(con_type[k] == PROJECTION)) {

               for (n = 0; n < Nptotal; n++) 
	            *(weightp + NN + n) = 1 ;
	  
               sprintf(message,
		    "\nReading Np target info for constraint # %d...", k+1) ;
               printTwice(message) ;

               readEpixFile(ta_filename[k], target + NN) ;

               if (strcmp (wt_filename[k], "full") != 0) 
	            readEpixFile(wt_filename[k], weightp + NN) ;

	       if (con_type[k] == TARGET) 
	            con_coef[k] *= volcoef * inv_mean_wt(weightp + NN) ;

               sprintf(message, "Initial analysis:") ; 
               printTwice(message) ;

	       analyzeTar(knownp, weightp + NN, Nptotal) ;
               NN += Nptotal ;
          }
       }
     }
}

real	inv_mean_wt(real *wtarray)
{
     int	n ;
     real	swsq = 0 ;

       /************************************************************
       We calculate  N / Sum(wt^2) ;
       ************************************************************/

     for (n = 0; n < Nptotal; n++, wtarray++)
	  swsq += *wtarray * *wtarray ;

     if (swsq == 0) 
	  Fatal("All weights are zero! - Quitting.") ;
     
     return (Nptotal / swsq) ; 
}

		/************************************************************
		Transfer the solution from this iteration (nump) to the 
		solution arrays. 
		************************************************************/
void updateKnownValues()	
{
     void	reviseMinMaxp() ;
     int        n ;
     real	sum ;	/* sum of snump	*/
     real	variance = 0 ;

       /************************************************************
       Add solutions from this iteration to solution set and 
       subtract same amounts from the ccg limiting arrays.
       ************************************************************/
     reviseMinMaxp(Npextended, nump) ;

     for (n = 0; n < Npextended; n++) {
     	*(snump+n) += *(nump+n) ;
        variance += fabs(*(nump+n)) ;
     }
       /************************************************************
       Provide diagnostics                              	  
       ************************************************************/
      sum = 0 ;
          
     for (n = 0; n < Npextended; n++)  
          sum += *(snump + n) ;
          
     sprintf(message,
       "Cumulative sum of recovered electrons is %g", sum) ;
     printTwice(message) ;

     variance = sqrt(variance/Npextended) ;

     sprintf(message,
       "Variance of electrons/voxel for this iteration is %g\n", variance) ;
     printTwice(message) ;
}
		/************************************************************
		Reinitialize nump; used to be in updateKnownValues(), but I
		found a need for nump between that call and the next outer
		iteration.
		************************************************************/
void reinitNump()	
{
     int        n ;

     for (n = 0; n < Npextended; n++) 
	*(nump +n) = 0 ;
}

		/************************************************************
	        Write out the current solution        
		************************************************************/
void writeSol(char *name) 
{
     void	analyzEpixFile(real *) ;
     void	writeEpixFile(char *, real *) ;
     void       write_HDFnump(char *, real *) ;
     void       write_HDFmap(char *, real *, real, real) ;
     void       transform() ;
     int        get_min_max(real *, int, real *, real *) ;
     void       printXmap() ;

     real       sum = 0 ;		/* sum of model + recovered electrons */
     int	k, n, NN ;

     /************************************************************
     The solution consists of the cumulative electrons found. 
     The known starting model is added in.  
     ************************************************************/

     for (n = 0; n < Npextended; n++) 
          *(totp+n) = *(snump+n) + *(knownp+n) ;

     for (n = 0; n < Npextended; n++) 
          sum += fabs(*totp + n) ;
          
     if ((strcmp (md_filename, "empty") != 0)) 
        sprintf(message,
	  "\nTotal electrons (starting model plus recovered): %g\n", sum) ;
     else
        sprintf(message,
	  "\nTotal electrons recovered: %g\n", sum) ;

     printTwice(message) ;
     
     printTwice("\t Analysis of electron densities:") ;

     analyzEpixFile(totp) ;

     for (k = 0, NN = 0; k < Nconstraints; k++, NN += Nptotal)  {

          if ((con_type[k] == TARGET)||(con_type[k] == PROJECTION)) {

               sprintf(message, "\nAnalysis over target # %d:", k+1) ; /* 1-origin */
               printTwice(message) ;
	       analyzeTar(totp, weightp + NN, Nptotal) ;
          }
     }
     writeEpixFile(name, totp) ; 

  /*******************************************************************
  Transform (in-place) .bin info to map.  Note: in addition to the 
  regular map (in xplor or hdf format), we also write an ordinary .bin
  file for other post-processors.  Note too that .bin files do NOT use 
  HDF format, regardless of user's input re. HDFOUT.
  *******************************************************************/
     
     transform(Nptotal, totp, 0) ;
     get_min_max(totp, Nptotal, &ar_min, &ar_max) ;
  
     if (HDFOUT) {
          write_HDFmap(name, totp, ar_min, ar_max) ;
     }
     else {
          sprintf(message," %7d %7d %7d %7d %7d %7d %7d %7d %7d",
             Nx, 0, Nx-1, Ny, 0, Ny-1, Nz, 0, Nz-1) ;
          printXmap(name, message, Nx, Ny, Nz, totp) ;
     }
     sprintf(message, "%s_map", name) ;
     writeEpixFile(message, totp) ;
}

void	analyzEpixFile(real *array)
{
     int	get_min_max(real *, int, real *, real *) ;
     int	n, j, jbunch, sj ;
     int	Ncount, nmax ;
     real	ar_min, ar_max ;
     double	decade[10], sdecade[10] ;	/*** 16/2/2005 NOTE double ***/
     real	value ;
     real	from, to, jfrom, jto ;
     int	nx, ny, nz ;

     for (n = 0; n < Nptotal; n++) 
	  if (fabs(*(array+n)) < ROUNDOFF)	/*CHANGED May 3, 2005 from
						1.e-10 to 1.e-5 */
	       *(array+n) = 0 ;

     nmax = get_min_max(array, Nptotal, &ar_min, &ar_max) ;

     /* Print message indicating range of solutions and position of maximum */
     sprintf(message,
             "Range is (%g, %g) el/voxel, (%g, %g) el/cub %s.\n",
             ar_min, ar_max, ar_min/incube_to_vox, ar_max/incube_to_vox,
             units) ;
     printTwice(message) ;

     nx = nmax % Nx ;
     ny = (nmax/Nx) % Ny ;
     nz = nmax / (Nx*Ny) ;
     if (nz >= Nz)
        nz -= Nz;		/* use simple lattice coords */

     sprintf(message, "Maximum  is at lattice point (%d %d %d)", nx, ny, nz) ;
     printTwice(message) ;

   /* Now for distribution of data ... */

     for (j = 0; j < 10; j++)
	  decade[j] = 0 ;

     for (n = 0, Ncount = 0; n < Nptotal; n++) {
          if( ar_max != 0 )
		  value = *(array+n)/ ar_max ;
	  else
	      value = 0.;

          j = (int) (value * 10.) ;
          if (j > 9) j = 9 ;
          if (j < 0) j = 0 ;
          decade[j] += 1 ;
          Ncount++ ;
     }

     for (j = 0; j < 10; j++)
	  decade[j]  *= 100. / Ncount ;

     /* Determine degree of bunching */

     for (jbunch = 10, j = 0; j < 10; j++)
          if (decade[j] > 50.)
               jbunch = j ;

     /* Print message indicating distribution of unbunched data */

     if (jbunch == 10) {
          printTwice
          ("Fractional distribution of electron densities") ;
          sprintf(message, "\t(based on %d out of %d entries):\n", Ncount, Nptotal) ;
          printTwice(message) ;
          printTwice("     Range         % \n") ;
          sprintf(message, "    < 0.1:     %#6.2f", decade[0]) ;
          printTwice(message) ;

          from = 0.1 ;
          for (j = 1; j < 9; j++) {
               to = (j+1) / 10. ;
               sprintf(message, "%3.1f - %3.1f:     %#6.2f",
                    from, to, decade[j]) ;
               printTwice(message) ;
               from = to ;
          }
          sprintf(message, "    > 0.9:     %#6.2f", decade[9]) ;
          printTwice(message) ;
     }

     /* Handle data with bunching.  */
     else {
          for (j = 0; j < 10; j++)
               sdecade[j] = 0 ;

          for (n = 0, Ncount = 0; n < Nptotal; n++) {
               if( ar_max != 0 )
		  value = *(array+n)/ ar_max ;
	       else
	          value = 0.;

               j = (int) (value * 10.) ;
               if (j > 9) j = 9 ;
               if (j < 0) j = 0 ;
               Ncount++ ;

               if (j == jbunch) {
                    sj = (int) (value*100 - 10*jbunch) ;
                    if (sj > 9) sj = 9 ;
                    if (sj < 0) sj = 0 ;
                    sdecade[sj] += 1 ;
               }
          }
          for (sj = 0; sj < 10; sj++)
               sdecade[sj]  *= 100. / Ncount ;

          printTwice(
          "Fractional distribution of electron densities > 0.5% of maximum") ;
          sprintf(message, "\t(based on %d out of %d entries):\n", Ncount, Nptotal) ;
          printTwice(message) ;
          printTwice("Full Range         %       Subrange            % \n") ;

          jfrom = jbunch/10. ;
          jto = jfrom + 0.01 ;
          sprintf(message,
            "    < 0.1:     %#6.2f     %4.2f - %4.2f:     %#6.2f",
                  decade[0], jfrom, jto, sdecade[0]) ;
          printTwice(message) ;

          from = 1 / 10. ;
          jfrom = jto;
          for (j = 1; j < 9; j++) {
               to =   from + 0.1  ;
               jto = jfrom + 0.01 ;
               sprintf(message,
             "%3.1f - %3.1f:     %#6.2f     %4.2f - %4.2f:     %#6.2f",
                    from, to, decade[j], jfrom, jto, sdecade[j]) ;
               printTwice(message) ;
               from = to ;
               jfrom = jto ;
          }

          sprintf(message,
                "    > 0.9:     %#6.2f     %4.2f - %4.2f:     %#6.2f",
                  decade[9], jfrom, jfrom + 0.01, sdecade[9]) ;
          printTwice(message) ;
     }
}

		/************************************************************
                Read in / write out binary file.
		************************************************************/

void	readEpixFile(char *name, real *array) 
{
     FILE 	*fp;
     float	*f_array ; 	/* for reading float */
     int	n ;
     int	Nread ;
     char	*ext = BINSUF ;
     void	read_HDFnump(char *, float *) ;

     for (n = 0; n < Nptotal; n++) 
          *(array + n) = 0 ;

     if (strcmp (name, "empty") == 0) {
         printTwice
     ("Note: 'empty' is a keyword signifying that all elements in the array are 0") ;
	  return ; 
     }

     if (strcmp (name, "full") == 0) {
         printTwice
     ("Note: 'full' is a keyword signifying that all elements in the array are 1.") ;

          for (n = 0; n < Nptotal; n++) 
               *(array + n) = 1. ;
          
	  return ; 
     }

     f_array = (real *) e_malloc(Nptotal*sizeof(real), "readEpixFile") ; 

     sprintf(solution_filename, "%s", name) ;

     if ((strstr(solution_filename, ext) == NULL))
	  strcat(solution_filename, ext) ;

     if ((fp = fopen(solution_filename, "rb")) == NULL) {
          sprintf(message, "readEpixFile: Cannot open %s", solution_filename); 
	  Fatal(message) ; 
     }
     Nread = fread(f_array, sizeof(float), Nptotal, fp); 
     if (Nread != Nptotal) {
          sprintf(message, 
     "Expected length : %d, actual length: %d", Nptotal, Nread) ;
          printTwice(message) ;
          sprintf(message, 
     "readEpixFile: wrong number of items in file %s", solution_filename); 
	  Fatal(message) ;
     }
        
     for (n = 0; n < Nptotal; n++) 
          *(array + n) = *(f_array + n) ;
     
     e_free(Nptotal*sizeof(float), f_array) ; 

     fclose(fp) ;
}
void	writeEpixFile(char *name, real *array) 
{
     void	write_HDFnump(char *, real *) ;
     FILE 	*fp;
     float	*f_array ;	/* for converting real to float  */
     int	n ;	
     char	*ext = BINSUF ;

     f_array = (float *) e_malloc(Nptotal*sizeof(float), "writeEpixFile") ;

     for (n = 0; n < Nptotal; n++) 
          *(f_array + n) = *(array + n) ;

     if (strstr(name, ext) == NULL)
          sprintf(solution_filename, "%s%s", name, ext) ;
     else
          sprintf(solution_filename, "%s", name) ;

     if ((fp = fopen(solution_filename, "wb")) == NULL)
     {
          sprintf(message, "writeEpixFile: Cannot open %s", solution_filename); 
	  Fatal(message) ;
     }
     if (fwrite(array, sizeof(real), Nptotal, fp) != Nptotal) {
          sprintf(message, "writeEpixFile: Error writing %s", solution_filename); 
	  Fatal(message) ;
     }
     fclose(fp);

     e_free(Nptotal*sizeof(float), f_array) ; 
}

void	analyzeTar(real *array, real *wt, int N)	
{
     int	n ;
     real	Sum = 0 ;
     real 	count = 0 ;
     char	messout[MAXSTRING] ;

     /* Find the average electron level (in el/voxel and el/cub units) over the 
     targetted region (for which *wt > 0).   */

     for (n = 0; n < N; n++) {
	  Sum += *(array+n) * *(wt+n) ;
	  if (fabs(*(wt+n)) > 0)
             count++ ;
     }

     sprintf(messout, 
     "\tAverage over targetted region = %7.3f el/voxel. %7.3f el/cub %s",
          Sum/count, Sum/(count*incube_to_vox), units) ;
     printTwice(messout) ;
}

int	get_min_max(real *array, int N, real *a_min, real *a_max)
{     
     int        n, nmax ; 

     *a_min = 10000. ;
     *a_max = -10000. ;
     nmax = 0 ;

     for (n = 0; n < N; n++) {
          if (*(array+n) < *a_min)
               *a_min = *(array+n) ;
          if (*(array+n) > *a_max) {
               *a_max = *(array+n) ;
               nmax = n ;
          }
     }
     return(nmax) ;
}
